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Abstract The problem of the generation of an intermediate 
image between two given images in an image sequence is 
considered. The problem is formulated as an optimal control 
problem governed by a transport equation. This approach 
bears similarities with the Horn & Schunck method for op- 
tical flow calculation but in fact the model is quite different. 
The images are modelled in BV and an analysis of solu- 
tions of transport equations with values in BV is included. 
Moreover, the existence of optimal controls is proven and 
necessary conditions are derived. Finally, two algorithms 
are given and numerical results are compared with existing 
methods. The new method is competitive with state-of-the- 
art methods and even outperforms several existing methods. 
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1 Introduction 

Image sequence interpolation is the generation of intermedi- 
ate images between two given images containing some rea- 
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sonable motion fields. It is mainly based on motion estima- 
tion and has broad applications in the area of video compres- 
sion. In video compression, the knowledge of motions helps 
remove the non-moving parts of images and compress video 
sequences with high compression rates. For example in the 
MPEG format, motion estimation is the most computation- 
ally expensive portion of the video encoder and normally 
solved by mesh-based matching techniques, e.g. blocking 
matching, gradient matching [37|. While decompressing a 
video intermediate images are generated by warping the im- 
age sequence with motion vectors. 

Another possibility of image interpolation is based on 
optical flow estimation. Since Horn and Schunck proposed 
the gradient-based method for optical flow estimation in their 
celebrated work |26|, this field has been widely developed 
till now. For example, instead of the linear constraint in the 
Horn & Schunck method one applies the non-linear isotropic 
constraints ll6l [T3l . anisotropic diffusion constraints 03011201 
and TV constraint l38l for preserving the flow edges, which 
is very useful for motion segmentation. Dealing with large 
displacements in image sequences one develops warping tech- 
nique lfT2l to estimate the flow field in a robust way. How- 
ever, in ll24ll is shown that the Horn & Schunck method is 
only suited for optical flow estimation, but not for matching 
image intensities, especially in case of large displacements, 
see also the argumentation in [34|. 

Borzf, Ito and Kunisch considered the optical flow prob- 
lem in the optimal control framework |10|. Due to an op- 
timal control formulation the estimated flow field is also 
suitable for image interpolation, since one searches the flow 
field such that the interpolated image has a best matching to 
a given image in the sense of some norm. In this paper we 
modify the model proposed in iflOll for interpolating inter- 
mediate images between two given images and analyze the 
well-posedness of the corresponding minimizing problem. 
In the end we introduce an efficient numerical method for 
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solving the optimality system and we also propose a modi- 
fication of the segregation loop of the optimality conditions 
system, which give better interpolation results and is robust 
with respect to the choice of regularization parameter. To 
evaluate our proposed interpolation methods we will utilize 
the image database generated by Middlebury College and 
compare our results using the evaluation method of Middle- 
bury with the results in P4ll . 

2 Modeling 

We are interested in finding a flow field, which is suitable 
for image matching. It means that instead of minimizing the 
optical flow constraint equation directly, we utilize the trans- 
port equation to fit a given image uq to another given image 
ut in the sense of some predefined norm in the cost func- 
tional. 

Let us model the optimal control problem governed by 
the transport equation. Consider the Cauchy problem for the 
transport equation in [0, T] x £2, Q. c R rf (generally d = 2): 

{d t u(t,x)+b(t,x) -V x u(t,x) =0 in]0,r] xQ, 
(1) 
u(0,x) — uq(x) in Q. 

Here b : [0, T] x Q. — > W 1 is an optical flow field, uq is a 
given initial condition and u is an unknown function depend- 
ing on t and x. We define the nonlinear solution operator of 

G.XxY — > Z, 
(u Q ,b) i ^ u, 

where X,Y,Z are normed spaces to be specified. Then, we 
define a linear "observation operator" Ej : u M- u(T), which 
observes the value of u at time T. By the chain (uo,b) n- 
u h-> u{T) we have the "control-to-state mapping" 

S:XxY — >U, 
S : (uQjb) i-> u(T). 

The space U is a subspace of Z, which not involves time 
t. The continuity of S will be investigated in the concrete 
contexts. Our intention is to find the flow field b such that 
the corresponding image S(uQ,b) matches the image uj at 
time T as well as possible. This motivates to minimize the 
functional A \\S(uQ,b) — ut\\u- However, this problem is ill- 
posed and an additional regularization term is needed. This 
regularized optimal control problem can be formulated as 
minimizing the following cost functional 

MJ(b) = h\S(uo,b)-u T \\l + ^\\b\\l, (2) 

ogF Z Z 

1 http : //vision. middlebury . edu/f low/data/ 



subject to divb = 0. (3) 

We use Tikhonov regularization to stabilize the cost func- 
tional and X is the regularization parameter. In the frame- 
work of optimal control [29 36| we call b the control and 
u the state. According to the conservation law lf25l and the 
divergence theorem [32], the divergence free constraint of 
b will make the flow volume conserving, smooth and vary 
not too much inside the flow field of a moving object. Such 
properties are desired to be enjoyed in image interpolation 
in case that the moving objects are not getting deformed. 
Such constraint is not new for optical flow estimation and 
was similarly introduced as a regularization constraint e.g. 

in (SEED. 

We emphasize, that our model is considerably different 
from the Horn & Schunck approach which is based on the 
optical flow constraint. There one has a given image u and 
a given derivative d t u (both at time fg) and one finds a flow 
field b = (v,w) by minimizing 

/ (d t u-b-Vu) 2 dt+ f \Vv\ 2 +\Vw\ 2 dx. 
Jn Ja 

The main conceptual difference between this approach and 
ours is that Horn & Schunck just consider one time to and 
match the flow field only to that time. Hence, it is unclear in 
what sense the produced field b could be useful to match a 
given image with another one. Our approach uses two given 
images and tries to find a flow field b which transports the 
first image as close as possible to the second image. The "op- 
tical flow constraint equation" now enters as a constraint to 
the optimization problem and not in the objective functional 
itself. 

In next chapter we will give some adequate spaces for 
u and b. Especially we are interested in images uq and uj 
which are of bounded variation. Hence, we introduce the so- 
lution theory of transport equations equipped with a smooth 
flow field and a BV image as initial value. Especially we 
need to work out conditions under which the BV -regularity 
is propagated by the flow field. Then, we will analyze the 
existence of a minimizer of problem (f2]i restricted to ([TJ and 

CD. 



3 Analysis of Well-posedness 

To analyze the solution operator G we use the method of 
characteristics. We start with the analysis of the correspond- 
ing ODEs, then derive existence results for initial values uq 
which are of bounded variation and finally derive a result on 
the weak sequential closedness of G. Together this shows 
the existence of an optimal control in the respective setting. 
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3.1 Basic Theory of ODE 

It is well-known that the solution theory of transport equa- 
tions has a tight relationship with the ordinary differential 
equation 

y(0 =b(t,y(t)) t€l, 

(4) 

7(a) = xq in £2. 

Regarding the solution theory of @, the existence and unique- 
ness of a solution can be derived by the theorem of Picard- 
Lindelof ll23l if b is Lipschitz continuous in space and uni- 
formly continuous in time. We can also relax the assump- 
tion on t of b to be integrable by the following Caratheodory 
theorem [4|, which a general version of the Picard-Lindelof 
theorem: 

Theorem 1 (Caratheodory) Define I = [a,c] and £2 is a 

bounded subset in M. d . Suppose b : I x £2 -> R d so that 

1. t —> b(t,x) is measurable in I for every x G £2; 

2. there exists C > with \b(t,x) —b(t,x')\ < C\x — x'\ for 
a.e. t E I and every x,x' G £2; 

3. b(t,x) — for a.e. t E I and every x G d£2; 

4. the function m(t) — \b(t,xo)\ is integrable in I for xq G 
£2. 

Then, there exists a unique solution "f : / — > £2 with 

i 

f(t)=x + J b(s,y*(s))ds 
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to the Cauchy problem (|4). 

As a consequence of the proof, the flow 7*(?) is absolutely 
continuous in [a,c]. Generally, if we consider the solution 
in [0, T] with T > c, we can restart y* at (c,7*(c)) until the 
unique continuous solution arrives at time T. The backward 
flow is the special case when the time t is smaller than the 
initial time a. 

Next, we want to choose an appropriate function space 
Y for b, which is suitable for the control problem. Accord- 
ing to [3 1 the space of Lipschitz functions is equivalent to 
W L °°(£2) d , if £2 is a bounded, convex, open set. According 
to lfT31l lower regularity of the flow field (i.e. b G W ' p with 
p < oo) does not preserve Z?V-regularity. However, the norm 
in W l °° is not well suited as a penalty term since it is dif- 
ficult to determine the necessary optimality conditions of b 
equipped with the L°°— norm. Thus, we assume additionally 
that the domain £2 enjoys the strong local Lipschitz condi- 
tion 1 and use the fact that Hq (£2) d is continuously embed- 
ded into W l '°°(£2) d under this assumption, when dim(X2) = 
2. Considering the divergence-free constraint on b we set 

H 3 Q dn (£2) 2 :={feHi(£2) 2 \dwf = 0}. 

Adjusting the assumption on the time of b in Theorem[T]and 
previous conditions on £2 we will assume that 



- £2 C K 2 is a bounded, convex, open set with the strong 
local Lipschitz condition 

- fcGL 2 ([0,r];// 3 ' div (f2) 2 ) 

throughout the paper unless otherwise stated. A proper choice 
for the space U will be discussed in Section [3~3l 

In order to formulate the solution of transport equation 
in a convenient way, we give the concept of classical flow 

oa. 

Definition 1 The classical flow of vector field b is a map 

<P(t,x) :[0,T]xQ — >£2 
which satisfies 



dt 



(t,x)=b(t,<P(t,x)) m]0,T]x£2, 



(5) 



, <P(0, x) —x in £2. 



A helpful property of <P will be given in the following corol- 
lary. 

Corollary 1 For every t G [0, T] the mapping <P(t, •) : £2 — > 
£2 is Lipschitz continuous and a diffeomorphism. 

Proof The injectivity can be derived from the uniqueness 
of the backward flow: If the flow <P starts from two points 
x\ 7^ X2 and arrives at some t at the same point <P(t,xi) = 
<P{t,xi) =x, the backward flow starting from (t,x) will be 
not unique. Regarding the surjectivity: for every point y G £2 
one can find a backward flow starting from (t,y) 



Y(t')=y + J b{s,y(s))d 



s = x G £2 , 



according to TheoremQ] In case t' — yields &(t ,x) =y. 

The Lipschitz regularity of <P is easily shown by the 
Gronwall's lemma. For details we refer to lfl6l . 

Since the Lipschitz continuity gives only the local C 1 - 
regularity, the ^-regularity of 4>(f,-) in £2 one can follow 
the results in 1161 . which states that if b has C 1 -regularity 
in space, then the flow <£(?,•) is also C 1 in space. In fact, 
Hq(£2) 2 is continuously embedded into C l (£2) 2 , and hence 
we derive the statement. □ 



3.2 Solution Theory of Transport Equations 

In this subsection we will consider the transport equation 
with the initial value uq in BV. The BV space is a natural 
space for images, since BV contains the functions with dis- 
continuities along hypersurfaces, i.e. edges of images |[3). 
However, the propagation of BV regularity is a delicate mat- 
ter. We formulate first the solution of transport equations 
with a smooth initial value: 
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Corollary 2 Let uq e C l (£2) and <P be a classical flow of 
vector field b. Then the transport equation ([TJ has unique 
solution 



u(t,x) = UQO<P 1 (*,•)(.*) 



(6) 



Proof Let us test (Q~|i along the characteristics denoted by 

(f,#(f,*)) 

= -^(t,0(t,x)) + b(t,0(t,x))-Vu(t,0(t,x)) 
at 

du . , .. d<Z> , . _ . ^ . . . 
= — (f , #(*,*)) + ■^■('.*) ' Vm (' , *(*>*)) 



(«(f,-)o0(f,x)). 



This implies that every solution is constant along the char- 
acteristics. Adjusting the initial value we derive © is a so- 
lution to (Q} and the uniqueness follows immediately from 
the uniqueness of flow <£>. □ 

Equipped with a non-differentiable initial value the classic 
solution (0 will not work. Next, we give the definition of 
the solution of transport equations in the weak sense. 

Definition 2 (Weak solution) If b and uq are summable 
functions and b is divergence free in space, then we say 
that a function u : [0, T] x £2 — >• K is a weak solution of 
(Q]) if the following identity holds for every function (p £ 

C-([0,T[xQ): 

T 

J J u(d,(p +b -V(p)dxdt = - J uo(x)(p(0,x)dx. (7) 
on a 

In Theorem|4]it will be shown that © is actually the unique 
weak solution of (Q} with uq G BV(£2). Before we are able 
to deal with the proof, we recall briefly the weak* topology 

of BV EKED, 



EV{Q) 



u :-^> u„ 



-> u and Du„ 



Du 



in) 



which possesses convenient compactness properties in the 
following theorem [3|. 

Theorem 2 Let (u„) C BV(£2). Then (u n ) converges weakly* 
to u in BV(£2) if and only if (u n ) is bounded in BV(£2) and 
converges to u in L l (£2 ). 

To prove that (|6j is a weak solution of (Q]) it is common to 
use the technique of mollifiers ll2D . In short, we smooth the 
initial value with a mollifier r) £ with variance e, let e con- 
verge to zero and investigate the convergence of the solution 
with a smooth initial value to a nonsmooth initial value. This 
will be done in next theorem. 

Theorem 3 Assume uq £ BV(£2),(p and (p~ are diffeomor- 
phisms and Lipschitz continuous in £2. Then, the sequence 
((uq * TJ £ ) o q>) converges to uqo q> in the weak* topology of 

w(n). 



Proof Let us verify first the L 1 -convergence of (uq *rj £ )oq> 
and set <p (x) = y 

\(uq * rj £ ) o (p(x) — uqo (p(x)\dx 



I 



= f\uo*T] e (y)-uo(y)\\det(Vcp l (y))\dy 
a 

< ||«0*T?e-"o|| i i( X 3) || det ( V <P _1 )|| L - (i 2)- 

Let L be the Lipschitz constant of (p^ 1 i.e. L= ||V<p~ Wvfa)*' 
then ||det(V<p _1 )|| LIX ,^ is bounded from above by 2L 2 . To- 
gether with the approximation property of mollifiers this 
gives the L l —convergence. Regarding the weak* convergence 
of Radon measures V(uq * 77 £ ) we observe that for every 
y/eq°(X2) 2 it holds 

V((uQ*T} e )o(p)\irdx 



(wo * r\ £ ) o (pdivy/dx 

(u Q *ri £ )(y)div(\ifo(p- l (y))\detV(p- l (y)\dy 

J J ri £ (y ~ s)u Q (s)dsdiv(\i/ o (p- 1 (y^detVy- 1 (y)\dy 
a n 

J J rj £ (y-s)dw(\j/o(p-' l (y))\detV(p- 1 (y)\dyu (s)ds 



a n 



rj £ * (div(y/o <p )|detV<p |) (s)uo(s)d, 



(8) 



Since <p~' is C 1 and Lipschitz continuous in £2, the con- 
volved term belongs to L 2 (£2 ) . Recall that in the two dimen- 
sional case BV(£2) is continuously embedded into L 2 (£2), 
then utilizing the approximate property of mollifiers implies 
that the equation ([8]) converges to 



n 



(*) 



div(y/o<p 1 (i))|detV<p (s)\u (s)ds 



div\ir(Z)u ((p(Z))d% 



n 



\j/D(u o <p) 



In (*) we applied the Gauss-Green formula for the BV func- 
tions ETj. □ 

Remark 1 Under the same assumptions of Theorem [3] one 
can derive from Theorem|2]that ((uq * T] £ ) o <p) is uniformly 
bounded in BV (£2) and converges to uq o (p in L l (£2), actu- 
ally also in LP (£2 ) with p < 2 due to the approximate prop- 
erty of mollifiers and the fact BV(£2) has a continuous em- 
bedding into L 2 (£2) in the two dimensional case. 
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Lemma 1 Assume that uq G EV(£2), <p(t,-) and <p~ l (t,-) 
are diffeomorphisms in £2 for every t G [0,T] and <p(-,x) is 
absolutely continuous in [0, T] for every x G £2. Define 

u e (t,x) = (u *ri £ )o(p(t,x). 

Then, u £ &C([0,T];BV(Q)). 

We skip the proof of LemmaQ] since it is a trivial result 
utilizing the substitution technique introduced in the proof 
of Theorem[3] Now, we are able to prove the existence and 
uniqueness of the weak solution of the transport equation 
(ID- 
Theorem 4 If uq G BV(£2), then there exits a unique weak 
solution 



u(t,x)=UQO<P (t,-)(x) 

o/O belonging to L°°([0,T];BV(£2)). 



(9) 



Proof Consider the transport equation with initial value uq 
convolved with mollifier rj £ 

d t u(t,x)+b(t,x) -V x u(t,x) = in]0,r] x £2 

u(0,x) — uq * r] £ (x) in £2. 

Corollary |2]implies that there exists a unique solution u e of 
the form 

u e (t,x) = (u *ri £ )o<p- l (t,-)(x). 
Let us define 
u(t,x) = UQ o <t>~ l (t, 

where u(t,-) G BV(£2) according to Theorem [3] for every 
t G [0, T]. Remark[T]gives that u £ (t , •) converges to u(t , ■) in 
L 2 (£2) and u £ (t , •) is uniformly bounded in BV(£2). And ac- 
cording to Lemma[T]this yields that u £ is uniformly bounded 
in L°°([0,T];BV(Q)), which is continuous embedded into 
L 2 ([0 7 T};L 2 (Qj). Hence, there exists a subsequence (u £k ) 
of (w £ ) such that 



u £k ^umL 2 ([Q,T]-L 2 (Q)) 



(10) 



and u G LT([Q, T];BV (£2))- Due to the weak convergence 
of k £ . in L 2 ([0,T];L 2 (£2)), one can derive for every <p G 
C?([0,T[xQ) it holds that 



The upper convergence is valid since b G L 2 ([0 ,T];L 2 (Q) 2 ) 
and thanks to ( fTOb . The lower convergence can be deduced 
from the property of approximate identity. The left equality 
is valid for a smooth initial value and smooth vector field. 
Hence, all of them imply the right equality. 

Regarding the uniqueness of weak solution it is shown in 
l2l that the continuity equation, which is equal to the trans- 
port equation in case drvb = 0, has a unique solution in the 
Cauchy-Lipschitz framework, i.e. b G L l ([0, T];W l '°°(M. d )). 
Definitely, it is also valid under our assumption of b. 

Because of the uniqueness of the weak solution the con- 
vergence of subsequence (m 6 . ) in the previous proof can be 
proceeded to the whole sequence (u £ ). □ 



3.3 Existence of a Minimizer 

The goal of this subsection is to complete the cost functional 
(|2|i with some reasonable norm and investigate the existence 
of a minimizer of problem (0. First of all, we give the norm 
of the penalty term of (0 w.r.t. b. According to [ 1 1 an equiv- 
alent norm of Hq is 



1/2 



\d a b\\h(ny> 



(11) 



We can easily find out that the seminorm (J n \ V Ab\ 2 dx) l l 2 
is actually another equivalent norm of Hq(£2) 2 , since it is 
equivalent to (fTTT i. For the regularity of b in time we can 
give the equivalent norm ofL 2 ([0,r];//3(X2) 2 ) 



T 

llL 2 ([0,r];// 3 (i2)2) =/ ll V4 H f .')l!i2(£>)4^- 
n 



(12) 



As discussed above, we assume that uq and uj are BV-functions. 
Hence, BV seems to be a proper choice for the space U. 
However, since BV is continuously embedded in L 2 (Q) for 
d = 2 we use U = L 2 (£2 ) (we discuss this choice in more 
detail in Section|4|. Hence, our cost functional is 



J{b)^2\\S{uQ,b)~u T \\ Ll{n) 



\\VAb(t,-)\\l 2{n)4 



dt. 



(13) 



/ / u £k [d,q> + b ■ Vq>]dxdt 
on 



- Ju *rj £k (p(Q,x)dx 
n 



f f u[d t (p + b-V(p]dxdt 
o h 



- J UQ(p(0,x)dx. 
n 



Lemma 2 If ((p n ) and ((p^ 1 ) ore sequences of diffeomor- 
phisms in £2 and the Jacobian determinant detV<p„ is uni- 
formly bounded in IT (£2) by the upper boundC. Then, ((uq* 
TJ £ ) o <p„) is uniformly bounded in BV(£2) w.r.t. n. 

Proof It is easy to check that (uq * rj £ ) is uniformly bounded 
in BV(£2) according to Theorem [2] and [3] Suppose that the 
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upper bound is C. Let us verify first the L 1 —norm by setting 

y = <Pn H*) 



\(uo*rie)°9„ l \dx 



\u *ri £ \\detV<p n (y)\dy 



<C J \u *rj e \dy 
n 

< CC||w || £ l( fl ). 

Regarding the variation norm by ||«o|| V ar(fi) := In \Duo\dx 
we have 



|V(mo*T7 £ )o<p„ l \dx 



|V(Mo*i?e)(y)||detV<p n (y)|rfy 
n 

<cj \V(u *r} e )(y)\dy 
a 

< CC\\uQ\\ var{n) . 



□ 



Lemma 3 If(b n ) is uniformly bounded in L 2 ([0, r];// 3 (X2) 2 ) 
and uq 6 BV(£2). Define u n£ = (mo * T] £ ) ° 't',^ 1 and u' n £ = 
u„ j£ (t). Then, there exists a subsequence {u nkt£ ) such that 
u n , e converges to some limit u £ in L 2 ([0,T];L p (Q)) with 
p < 2 and weakly to u £ with p = 2. u* £ converges to u £ (t) 
in L P (Q) with p < 2 and weakly to u £ (t) with p — 2. 

Proof Recall that for every b„ there is a corresponding <P n 
s.t. n (t, •) e W l ^{n) 2 and || V0 n (t, •) || L > (fl)4 = Lip(<£„(f, •))• 
The Lipschitz continuity implies via Gronwall's lemma 



Lip(<£„(? , •)) < exp / Lip(b n (s, -))ds 



(14) 



Now we consider the minimization problem 

inf J{b) (15) 

Z>GL 2 ([0,r];fl 3 ' dlv (r2)2) 

with / according to (TT31 . Proving the existence of minimiz- 
ers is usually achieved by the direct method [7 1 and the most 
difficult part lies in the weak sequential closeness of the so- 
lution operator G with respect to b. 

Theorem 5 (Weak sequential closeness) Suppose the se- 
quence (b n ) £ L 2 ([0, r];//Q ,dlv (X2) 2 ) is uniformly bounded 
and converges weakly to b in L 2 ([0, T];H 3 (Q) 2 ). Let u n be 
the corresponding weak solutions of ([TJ with flow field b„ 
and initial value uq ( i.e. u n — G(uo,b) ). Suppose that u„ con- 
verges to u in L 2 ([0,r];L 1 (X2)) and u 6 L 2 ([Q, T];L 2 (Q)), 
then u — G(uQ,b). 

Proof Since (b„) converges weakly tofoinL 2 ([0,r];// 3 (X2) 2 ), 
it is also valid that 



b„ ^ b in L 2 {[Q,T];L 2 (Q) 2 ). 



(16) 



Let us consider the difference u„ — u applying a test function 
<peC?([0,T[xQ): 



1 

J J u„(d t (p + b n V(p)-u(d t (p + bV(p)dxdt 
o n 

T T 

d,(p(u n — ujdxdt + J J V<p ■ (u„b„ — ub)dxdt 
on on 



The boundedness of (b n ) inL 2 ([0, T];H 3 (Q) 2 ) gives the up- 
per bound of (TBI . Hence, the Jacobian determinant det Vc£>„(f , •) 
is also uniformly bounded in L°° (Q ) . According to Lemma|2] 
this implies that u' n E is uniformly bounded in BV(Q) w.r.t. 
n. Then, there exists a subsequence (hL e ) of (u' n £ ) such 
that u* £ converges to u' £ in L P (Q) (weakly for p = 2) with 
p < 2. Considering the integral over time one has 



Part (z) converges to zero, since u„ — > u inL 2 ([0, T];L l (12)). 
Regarding part (h) we can derive 

T 

J j V(p(u„b„ — ub)dxdt 
n 

T T 

V(pb n (u„ — u)dxdt + J j V(pu(b„ — b)dxdt 
on on 

< H V( PllL°°([0,7-]xi2) 2 ll^«llL 2 ([0,7-];L°°(i2) 2 ) \\ u n ~ Hl 2 {\0,T\,L 1 {O)) 




1 

- J J V(pu(b„—b)dxdt 



o n 



lim 



*n k ,e 



\u>{n) 



dt 



lim 

n k -*oo 



*n k ,e 



\LP{n) 



with p < 2. The exchange of the limit is valid since the inte- 
grand is bounded and with the same argument one can derive 
the weak convergence of u„ k:£ in L 2 ([0, T];L 2 (Q)). □ 



Since (b n ) is uniformly bounded in L 2 ([0,T];H i (Q) 2 ), it is 
also uniformly bounded in L 2 ([0, T];L°°(Q,) 2 ). Due to the 
convergence of u„ in L 2 ([0,T};L l (£2)) and ( [ToT l imply the 
two summands of last inequality converge respectively to 
dt -> zero. 

Since («„) are weak solutions of (Q]), the limit u is also a 
weak solution of (Q]), i.e. u = G(uo,b), □ 



Theorem 6 (Existence of a minimizer) Suppose uq G BV (£2 ), 
then the minimization problem ( 1151 l has a solution. 
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Proof Let [b„) C L 2 {[0, T];H ' "(Q) 2 ) be a minimizing se- 
quence of the cost functional. The coercivity of (TT~3T > is a 
natural property subject to the norm (fT~2l >. From the coerciv- 
ity one has (b„) is uniformly bounded in L 2 ([0, r];// 3 (f2) 2 ), 
then there is a subsequence (b„ k ) of (&„) converging weakly 
to b in L 2 ([0, r];// 3 (X2) 2 ). For each b„ there exits a unique 
flow <P~ , which is a diffeomorphism in and absolutely 
continuous in [0,T]. Define 

"n,e = (mo*^) ^ -1 - 

According to Lemma [JJ there exists a subsequence (u„ ki e), 
which converges to u £ E L 2 ([Q, T];L 2 (Q)) inL 2 ([0, T];L l '{Q)) 
and converges for every t £ [0, T] weakly to u £ (t) in L 2 (Q). 
Theorem[5]implies that u £ = (uq * rj £ ) o <£> 1 . Hence, it yields 
that 



fu'(pdx 

Q 

'21 



J u' £ (pdx 
h 



f u'(pdx 
h 



for every <p £ L (Q). The left and right convergences in the 
diagram are valid due to the property of approximate identi- 
ties according and then u' = uq o <J> 1 (f , ■). Hence, wL con- 
verges weakly to u' in L 2 (Q) for every f £ [0, T]. 

The l.s.c. of the first term in Sl3i can be easily derived 
from — ut — 1 u T — ut in L 2 (i2). And the l.s.c. of the 
second term in ([131 ) is valid due to the norm-continuity of 
b. □ 



4 First-order Optimality Conditions System 



5 Algorithms 

In this section we will present an efficient numerical algo- 
rithm to discretize the optimality conditions system. Regard- 
ing the forward and backward transport equations in ([T8l 
one can take advantage of explicit formula © and estimate 
the backward flow by the fourth-order Runge-Kutta method. 
Another possibility for solving the transport equations is to 
utilize the explicit high-order TVD schemes with flux limiter 
"superbee" [25,28, 10 1. It works very well for preserving the 
edges of images and avoiding oscillations of solutions. The 
last equation of (TT81 is a triharmonic equation which stems 
from the use of space Hq as penalty term in (QjJ. There are 
little articles about its numerical schemes, e.g. 1171 . But the 
algorithms are either not efficient or difficult to be applied 
directly. The motivation for this term was that b has to be 
Lipschitz continuous to obtain a unique flow <P. If we ap- 
ply some smooth initial flow b° in the discrete form of ( [T81 
and replacing A 3 with A in ([T8l still leads to smooth enough 
b. Actually, according to [18] an initial value uq, GL 2 (Q) is 
transported into an L 2 (£2) -function by a flow field b £ H l . 
Hence, in our context we can also work with the optimality 
system 



u t +b-Vu = 0, u(0) = uq 

Pt +b-Vp = 0, p(T) = -(u{T)-u T ) 
divb = 0, 



(19) 



XAb + Vq = pVu 1 b = 0,V„b = 0, 
Ab = on d£2. 



We use the Lagrangian technique to compute the first-order 
optimality conditions of control problem < TT~3T > governed by 
(Q]l and ([JJ. Let us define first the minimizing functional with 
Lagrange multipliers (p,q) 



o a 



o n 



We remark that the assumption uq.uj G BV is not present 
in this model anymore. One could easily use U = BV and 
the BV-norm for the difference u(T) — uj since this would 
only affect the right hand side of the adjoint equation. How- 
ever, in this case we have to ensure that the flow field b is 
Lipschitz- continuous. In numerical experiments we found, 
T t that this did not alter the results too much and hence, we use 

L(u,b,p,q)=J(u,b)+ / / (u,+b-Vu)pdxdt+ / / divbqdxdt , the optimality system ([19}. 

' ' The hierarchical processing according to |9|,i.e. acoarse 

to fine calculation, provides a good choice of b°. The qual- 
ity of b° depends strongly on the downsampling and upsam- 
pling procedures of images. 

With a divergence free initial value b° we propose a seg- 
regation loop in the spirit of iflOl to interpolate the interme- 
diate image at time t : 
Segregation loop I. 

Suppose n = !,■■■ ,Ni oop and Ni oop is the iteration number. 
Given uq,ut, b n ~ it), A . The iteration process for solv- 
ing ( [T9l at iteration n proceeds as follows: 



(17) 



the variable p is the adjoint state of u and q is the adjoint 
state of b. The functional derivatives of (ITTb w.r.t. u and b 
yield the first-order necessary conditions system 



u(0) = uq 

p(T) = -(u(T)-u T ) 



u t + b ■ Vu = 0, 
p t + b-Vp = 0, 
divb = 0, 

XA 3 b + Vq = pVu, b = 0,V n b = 0, 
Ab^QondQ. 



(18) 



1. Compute u n - x (t)yu"- l (t) and u n ~ l {T) by the forward 
transport equation using uq and b n ~ l . 
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2. Compute p n 1 (?) by the backward transport equation us- 
ing -(u"- l (T)-u T ) and b"' 1 . 

3. Compute b n (t) by the Stokes equation with right-hand 
side p"- l {t)Vu n - l (t) and a A". 

After N i ol ,p iterations the intermediate image u Nlo °p(t) ap- 
proximating u at time t. Moreover, we use a monotonically 
decreasing sequence (A"), which converges to a final A*. 
However, thanks to the theory of Stokes equations 11221 . we 
know that 



ll*(OII«i(fl) < X H/»(0V«(0llfl-i(fl) , a.e. t e [0,T]. (20) 



In practice we find out that if we choose (A") such that the 
norm of the right-hand side of (120b is monotonically increas- 
ing, the value of b(t) will be also increasing. However, the 
final A* cannot be chosen too small such that the minimizing 
process of (13[ is ill-posed. 

Moreover, since the system ( TT8l > is a necessary condi- 
tion of minimizing functional d X 3T >. one expects that the term 
\\u(T) — Ut\\i2/q\ is not very small. But since this is one of 
our final goals, we propose a modification of segregation 
loop I, which poses no requirement for choosing a specific 
sequence (A") and gives better approximation of intermedi- 
ate images. We modify segregation loop I as follows: 
Segregation loop II. 

Suppose n = !,■■■ ,Ni oop and Ni oop is the iteration number. 
Given uq,ut, b n ~ l (f ) , A . The iteration process at iteration n 
proceeds as follows: 

1. Compute u n - 1 (t),Vu n - l (t) and u n - l {T)by the forward 
transport equation using uq and b . 

2. Compute p"~ l {t) by the backward transport equation us- 
ing -(u"- l {T)-u T ) and b"' 1 . 

3. Compute the solution of the Stokes equations with right- 
hand side p n_1 (f)Vn* -1 (f) and A. Then, denote it by 
8b"- 1 {t). 

4. b"(t)=b n - 1 {t) + 8b"- 1 (t). 

In segregation loop II we utilize the system ([T9T i to estimate 
the update of the flow field and update the flow field in step 
4. This point of view is different from the original problem 
$1% , but interestingly this modification actually solves the 
necessary condition of another minimizing problem. If the 
segregation loop II converges, then the update 8b"- 1 (t) con- 
verges to zero. Since the initial value b° is divergence free 
and in each iteration the update flow 8b"- 1 is divergence 
free, the limit of b" is also divergence free. 

We denote u* ,p* ,b* ,q* the limits of particular sequences 
and in this case 8b* = 0. Setting the limits into (T% we de- 



( u* + b* -Vu* = u*(0)=u 

P ; + b*-Vp*=0 p*(T) = -[u*[T)-u T ) 

divb* =0 b* = on dQ 

[ V<?* = p*Vu* 



(21) 



Actually, (l2TT> is the optimality system of another constrained 
minimization problem, namely 



-IKCO-wrlltf^) 
subject to 

u*+b*Vu*=0 u*(0)=u 
div6*=0 b*=0ond£2. 



(22) 



(23) 



Compared to (fT3l the functional (l22l is not regularized. But 
if we stop the segregation loop II on time, i.e. the interpola- 
tion error does not vary too much, then it is not surprising 
that segregation loop II gives good approximation results of 
intermediate images. From the point of view of regulariza- 
tion theory, one may see the segregation loop II as a kind 
of a Landweber method for minimizing \\u(T) — ht||^/q\ 
which is inspired by a Tikhonov-functional. 

In the most cases the forward interpolation from uo to 
uj and the backward interpolation from uj to uq are com- 
plementary, since the flow is only able to transport objects 
from somewhere to somewhere, but not able to create some 
new objects. If in the forward case some new objects ap- 
pear, then in the backward case the new objects disappear. 
It means that backward interpolation is more suitable for in- 
terpolating the intermediate images. In practice, we take the 
average of forward and backward interpolations. 



5.1 Hierarchical Method 

In order to get a start value b° for the optimality system, 
the hierarchical processing is a good ansatz. It can be under- 
stood in level I in the following steps: 

1. Downsample the images into level /. 

2. Solve system ( fT9b in level / out and get b l . 

3. Upsample the optical flow into level I — 1 and get b'^ 1 . 

The estimated optical flow b l ~ l is a start value of the hier- 
archical method in level I — I. In coarsest level we assume 
the start value is zero. As above mentioned, the down- and 
up-sampling methods are decisively, i.e. it is supposed to 
lose the local structures of objects as small as possible while 
down- and up-sampling the images or the optical flow. 
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In practice, we apply bicubic interpolation [31] for the 
sampling, since it has fewer interpolation artifacts than bi- 
linear interpolation or nearest-neighbor interpolation. Com- 
pared to the Gaussian pyramid lfl4l the downsampled im- 
ages by bicubic interpolation does look not so blurred. 



5.2 Numerical Schemes for Transport Equations 



Another way for solving the transport equation is to utilize 
the characteristic solution. From © we know the keypoint 
is to solve the backward flow starting from (t,x) 



b(s,&) m[0,t[xn, 




(24) 



To discretize the transport equations we can use the second- 
order TVD scheme. It is also suitable for the backward trans- 
port equation, since we can reform it into the forward prob- 
lem by setting t' := T — t: 

Pf -b-Vp = 0, p(0) = -(w(0) -Mr). 

Suppose the image size is N x M, h and At are the mesh 
sizes in space and time, respectively with mesh index i = 
I,'" ,N,j — 1, • • • ,M in space and k— 1, • • • ,K in time. The 
stability condition of the scheme, usually called CFL condi- 
tion 1 3, is 

At 

Ocfl :=max(|v| max , \w\ max )— < 1. 

n 

by setting b := (v, w). In practice we choose At such that 
Ocfl = 0.1. The TVD scheme of the forward transport equa- 
tion is: 



U t \ 



VUr 



At 



•J 

h 



:(^3 ,.)' 



1 x[r 3 ,< 



•j 
h 

k 

i+ 1 J ' 



1 l*( r ;+ 3 P 
Z '+iJ 2 r 



where = max(vy,0), v ; ■ = min(v, ; -,0) and the flux differ- 
ence ratios are defined as 



u i+hi~ u ij 



u 



-l, 



■4 



"2J 



r - 3 ■ 



In the similar way we can discretize the term — wu y . The 
superbee limiter function is given by 

X(r) = max(0,min(2r, l),min(r,2)). 

To compute the spatial derivatives of images we use the stan- 
dard three-point formula: 

pu x \ij = 2f.{-Ui-ij + u i+ \j)pij, 
1 , 

PUy\ij = —{-Uij-l+U iJ+ i)pij. 



To solve d24] > numerical efficiently we use Runge-Kutta 4th 
order method PP . We discretize [f ,0] with time step At = 
0.1 and utilize a constant flow b over [t ,0] due to saving the 
memory and computational cost. In this scheme we have to 
interpolate the flow b(t,x) with some non-integer x, since 
only the flow b(t, •) with integer coordinates is given. For 
this we use bilinear interpolation (a bicubic interpolation 
leads to almost the same results). Then, we warp the im- 
age mo with the coordinates calculated by (124-b using cubic 
spline predefined in Matlab to approximate u(t,x). 



5.3 Finite Element Methods for Stokes Equations 

As previously mentioned, after replacing A 3 with A it is im- 
mediately seen that the last two equations in ( fT9l ) are the 
Stokes equations. Stokes flow estimation was investigated in 
ll33l and Suter applied the mixed finite element method [35] 
for solving it. Moreover, the approximation of velocity field 
b(t,-) and pressure q(t,-) will achieved by the polynomial 
of second order (P2) and first order (PI), so-called Taylor 
and Hood elements [ 19 1. If the chosen finite element spaces 
satisfy the inf-sup condition, also called LBB condition lfl9l 
[TTIl . then the method is stable. 

The variational problem of the Stokes equations reads as 
follows: 



a(b(t),v)+c(v,q(t)) = (f(t),v) Vv e V, 
c(b(t),w) =0, VweW 

and the bilinear forms are defined by 

a(b(t),v) = / XVb{t)Vvdxdy, 



(25) 



c(v,q(t)) = J (divv)q(t)dxdy, 
h 

(/(?), v) = -Jf(t)vdxdy, 
n 

where /:= pVu,V :—Hq(Q) 2 and 



W := < weL 2 (Q) 



wdxdy = 



10 



The discretization of (l25T > using the mixed finite element 
produces a linear system of the form 



co)\pq)~\ 



(26) 



The approximation coefficients bMN,PQ and /mn are w.r.t. 
the basis of finite element spaces V/, and Wft. The stiffness 
matrix A has the following block form: 



Ai 
Aj 



where A i = ( J n V(p i V(pjdxdy) j . , i, j = 1 , • • • , MN and <p, are 
the basic functions of V/,. The matrix C r has also a block 
form 



a 



Cj = { / ^-xjfjdxdy 



i = l, 



,MAr;; = l,...,(2 



Xj/jdxdy 



i = l,---,MN;j = l,---,Q 



Similarly, y/, are the basic functions of W/,. The vector / = 
{fiifif is composed of scalar products (/i,<Pi) and (/2,<Pi) 
for / = 1, • • • ,MAf. We derive the interpolation polynomial 
of /i,/2 w.r.t. the basic functions 



MN 



/2 = E/2(^)<P«, 



;=1 

where x, is the corresponding measurement point of <p, . Then, 

MAT » 

/i = ifi,9i) = L/i (■*/) / (PjVidxdy, i = 1, • ■ ■ ,MJV 

;'=1 ^ 

,. 

/i = (/I , <P«) = E /2 (*;) / ¥!/ WMy. / = AflV + 1 , ■ ■ • , 2MAT. 

For simplifying the estimation we just need to define the 
basic functions of a single element, i.e. a triangle or square, 
and derive the corresponding element stiffness matrix and 
element mass matrix, then assemble them into A\, C\, C2, 
/mn- 

Since the matrix in (l26l i is sparse and symmetric, but not 
positive definite, the system (l26l i can be numerically solved 
by the routine bicgstab predefined in MATLAB. 



6 Numerical Experiments 

6.1 Parameter Choice Rule 

The essential parameters of the quality of image interpola- 
tion are the regularization parameter A and the downsam- 
pling level /. Experimentally, we find out that the optimal 
regularization parameter X opt and I are coupled. The down- 
sampling level should be so adapted that at the lowest level 
L the estimated optical flow is accurate with a X^ pt . At the 
higher level I with I < N the parameter k l opt is larger than 
Xg pt . In practice, we choose X opt with I < N by the follow- 
ing strategy: 

1 . Find a pair (X^ pt , L) experimentally at the lowest level L. 

2. Choose l ( '-J such that X\ w } /X' opt £ [10 2 10 5 ] and the 
interpolation errors decrease at level I — I. 

The difference between segregation loop I and II lies in that 
segregation loop II equips with a constant X l opt at each level 
and segregation loop I applies a monotonically decreasing 
sequence converging to X l opt at each level. In case the image 
size is around 600 x 400 we set the lowest level L = 3 and 

A^e[io 5 io 5 - 5 ]. 



6.2 Numerical Results 

To illustrate the effect of our intermediate interpolated im- 
ages, we apply the interpolation error (IE) introduced by 1 8 1 . 
Moreover, the IE measures the root-mean-square (RMS) dif- 
ference between the ground-truth image u and the interpo- 
lated image u 

1 

N M 



where M xN is the image size. We test our methods on the 
datasets generated by Middlebury with public ground-truth 
interpolation: 

- Dimetrodon with size 584 x 388 

- Venus with size 420 x 380 

Every dataset is composed of three images and the mid- 
image is the ground-truth interpolation at time 0.5 if we as- 
sume the evolution process of three images lasts time T = 1 . 
To evaluate the interpolation we can compare our interpo- 
lation results with the ground-truth by means of IE mea- 
sure. The ranking of the interpolation results calculated by 
segregation loop I and II refers to Table Q] As in [ 8 1 men- 
tioned the Pyramid LK method and Mediaplayer™ are sig- 
nificantly better for interpolation than for ground-truth mo- 
tion, since e.g. Mediaplayer™ tends to overly extend the 
flow into textureless regions, which are not significantly af- 
fected by image interpolation. According to Table Q] segre- 
gation loop II works better than some classic methods and 



11 



more accurate than segregation loop I. The places where the 
interpolation errors take place refer to Fig. Q]-[2] As a result 
our methods especially segregation loop II work efficiently 
in image interpolation. 





Dimetrodon 


Venus 


Segregation loop I 


2.25 


6.67 


Segregation loop II 


1.95 


3.63 


Stich et al. 


1.78 


2.88 


Pyramid LK 


2.49 


3.67 


Bruhn et al. 


2.59 


3.73 


Black and Anandan 


2.56 


3.93 


Mediaplayer™ 


2.68 


4.54 


Zitnick et al. 


3.06 


5.33 



Table 1 Interpolation errors calculated by our methods using the Mid- 
dlebury datasets by comparison to the ground truth interpolation with 
results taken from 1 34 1 . 

The whole interpolation process of Middlebury datasets 
is accomplished by 9 generated images respectively using 
segregation loop I and II. The additional data generated into 
films are given in Online Resource. 

7 Conclusion and Outlooking 

The approach to image sequence interpolation by optimal 
control of a transport equation has proven to be useful and 
competitive to existing methods. While we started to model 
the images in BV we ended up with an algorithm which does 
not exploit this regularity but merely uses the L 2 -structure. 
This was due to the fact that one needs Lipschitz-continuous 
flow fields to preserve BV -regularity lfT31 . Hence, we finally 
used H l flow fields. However, this still imposes some reg- 
ularity on the flow field and discontinuous flow fields are 
still not allowed. In further work it may be interesting to use 
BV vector fields and hence try to transport an image with a 
possibly discontinuous flow field. Another open question is, 
how to deal with objects which appear in the second image 
but are not present in the first image. One possibility could 
be to use heuristic techniques to estimate motions which oc- 
clude or disclose objects as described in [34|. 
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Fig. 1 (a) uq. (b) uj. (c) «o plus the colored difference between uq and uj. (d) The groundtruth interpolation at time T/2 from the Middlebury 
datasets. (e) The generated interpolation at time T/2 by segregation loop I. (f) The absolute difference between (d) and (e). (g) The generated 
interpolation at time T/2 by segregation loop II. (h) The absolute difference between (d) and (g). 
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Fig. 2 (a) Uq. (b) uj. (c) uq plus the colored difference between «o and uj. (d) The groundtruth interpolation at time T /2 from the Middlebury 
datasets. (e) The generated interpolation at time T /2 by segregation loop I, (f) The absolute difference between (d) and (e). (g) The generated 
interpolation at time 7/2 by segregation loop II. (h) The absolute difference between (d) and (g). 
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